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SATELLITE  BREAKUPS:  SURVEY  OF  A DYNAMICAL  THEORY 
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INTRODUCTION 

The  phenomenon  of  artificial  satellite  breakups,  or  sudden  dis- 
integrations, presents  a plethora  of  interesting  dynamical  problems 
concerning  the  motion  of  the  fragments.  The  basic  problem,  which 
will  be  called  the  direct  problem,  is  to  predict  the  motion  of  the 
ensemble  of  particles  from  knowledge  of  the  characteristics  of  the 
disrupting  mechanism,  specifically  the  relative  velocity  imparted  to 
the  fragments.  An  important  variation  on  the  basic  direct  problem  is 
the  continuous  source  problem  which  considers  particles  dispersing 
over  a finite  interval  of  time  rather  than  instantaneously.  The 
eventual  distribution  of  the  particles  after  a long  period  of  time, 
the  so-called  asymptotic  problem,  is  also  of  interest.  Finally,  one 
may  invert  the  problem  and  seek  to  determine  the  characteristics  of 
the  disrupting  mechanism  from  observations  of  the  particle  distribution 
at  some  time  after  breakup  — the  so-called  inverse  problem. 

Each  of  the  problems  Just  mentioned  has  its  place  in  the  study  of 

satellite  breakups.  However,  their  applicability  does  not  end  there 

because  several  astronomical  phenomena  may  be  mentioned  which  are 

represented  by  such  a dynamical  system.  For  example,  there  are  the 

problems  of  meteor  streams  originating  from  a disintegrating  comet 

(Jacchia,  1963),  expanding  stellar  associations  (Blaauw,  19?2)  and 
Note:  Manuscript  submitted  June  7,  1978. 
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fragmented  asteroids  (Wiesel,  1977). 
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I 

The  purpose  of  this  report  is  to  survey  the  dynamical  theory  of 
satellite  breakups  as  it  has  been  developed  over  the  past  two  years 
or  so  (Heard,  1976;  1Q77  a,b).  In  view  of  the  variety  of  possible 
applications  of  the  theory,  each  aspect  of  it  has  been  developed  in 
as  general  a form  as  possible.  This  is  not  intended  to  distract  the 
development  from  its  primary  motivation,  namely  satellite  breakups, 
but  rather  to  make  it  accessible  to  the  largest  possible  sphere  of 
application.  Ibis  report  consists  of  some  previously  published 
results,  some  generalizations  of  previously  published  results,  and 
some  new  material.  The  goal  is  to  give  a complete  and  coherent 
presentation  of  the  theory  as  it  now  stands. 

The  basic  method  underlying  each  facet  of  the  theory  is  to 
represent  the  finite  number  of  discrete  fragments  as  a continuum  of 
particles  distributed  in  phase-space.  This  approach  is  borrowed  from 
the  field  of  stellar  dynamics,  where  it  is  known  as  the  continuum 
theory  ( Contopoulos , 1066).  In  this  approach,  all  results  are  derived 
from  solutions  of  the  phase-space  conservation  equation  ( Chandrasekhar' , 

I960). 

DIRECT  PROBLEM  (instantaneous  Disruption) 

Mathematical  Formulation  and  Formal  Solution 
Consider  an  ensemble  of  non-interacting  particles  moving  in  a 
common  force  field  such  that  the  equations  of  motion  for  an  individual 
particle  are 


I t 


II 


•j  “ v i a.  5 , * ) 


s = - & . " 


where  the  n-vectors  q ,2  are  the  coordinates  and  momenta,  respective- 
ly. Let  be  the  phase  space  density  function  for  the 

ensemble  such  that  a volume  at  point  ^ ,£  contains  dN 

particles  at  time  t where 

According  to  Chandrasekhar  (i960),  the  distribution  function  f 
satisfies  the  following  first-order,  linear  partial  differential 


equation 


57  0 


if  there  is  no  source  of  particles  present.  Equation  (2)  may  be 


written 


2i  = - f A 

Dt 


where  JL  denotes  the  ’’Stokes  derivative" 
Dt 


-L  . _L  - Y • V * Y -V 

Dt  ~ at  - % - t 


♦ ^1 

fr;Ln. 


Introduce  propagator  functions 


9-rH>  (1*) 

Defined  such  that  a particle  at  at  time  t will  be  at  3rt^,£.t) 

5 att ime  t * r . 

The  formal  solution  of  (3)  which  satisfies  the  initial,  condition 


^*,0)  = 


may  now  be  written 


■ « 


where 


J it  ; 


and  it  is  understood  that  the  integral  (7)  is  to  be  evaluated  along 
the  trajectory  passing  through  \ ^ at  time  t . 

A fundamental  descriptor  of  the  ensemble  is  the  spatial  density 
t)  . It  is  obtained  by  integrating  r*  over  all  momenta,  i.e.. 


f * J • 


The  spatial  density  specifies  the  apparent  shape  and  dispersion  of  the 
particle  cloud  at  any  time.  The  quadrature  (8)  is  easily  performed 


h 


-ormaily  it  the  particles  emanate  from  a common  point. 

To  describe  dispersion  from  a common  point,  j , we  may  use  the 

7 * 

following  initial  condition: 

F(V  J , (9) 


where  5 ( ) denotes  the  Dirac  delta-function  and  G-  is  an  arbitrary 

function  which  specifies  the  initial  distribution  of  momenta. 

The  integral  (8)  is  evaluated  by  applying  two  theorems  from  the 
calculus  of  S-  functions  (Friedman,  1956).  Theorem  I relates  a 
multi-dimensional  & -function  to  the  product  of  one-dimensional 
8 -functions,  and  Theorem  II  specifies  the  behaviour  of  S -functions 
under  a change  of  variable: 


£-£.}=  s(xi- ... 


11  * J ru)  dx  _L  f(j  j ; 

where  is  determined  from  the  equations 

K <*,)  = ° , i a j 


and  J is  the  Jacobian  determinant 

J= 

3(»,,  ■ • 


evaluated  at  x = x 
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Now,  substitute  (9)  and  (6)  into  (8)  to  obtain 


\\ 


H 


According  to  Theorems  I and  II  we  then  obtain 

* Ijj  &[£-t(V**..o]«f{-rtv>irlo}  J (11) 

where  ^ is  the  solution  of  the  equation 

= %*  j (12) 

and 

T = g-t,i  ■ •••  j (13) 

3 ( P,  > '••jPJ 

Thus,  the  calculation  of  the  evolution  of  spatial  density  is 
reduced  to  the  evaluation  of  a Jacobian  determinant  (13)  and  the 
solution  of  an  equation  (12)  which  is  usually  transcendental. 

Solution  for  Small  Velocity  Increments 

In  virtually  all  of  the  applications  envisaged  for  the  theory, 
the  relative  velocities  of  the  particles  are  small  in  comparison 

with  the  total  velocity  of  the  parent.  The  mathematical  import  of 

this  is  that  the  equations  of  motion  may  usually  be  linearized. 

With  this  motivation,  we  now  reduce  the  general  results  of  the 


A 


previous  section  to  the  special  form  assumed  for  linear  systems. 


J-l 


Consider  a dynamical  system  for  which  the  equations  of  motion 
of  an  individual  particle  are 


i sA  x ♦ f 


(lU) 


where  A is  a uou-singuiar  matrix  and  f-  is  a vector.  Let  Y(t)  be 
a matrix  solution  of 


IT  * A* 


(15) 


whose  columns  are  linearly  independent,  and  let 


lit.O  =Y(t)  . 


(16) 


The  matrix  $ is  the  matrizant  of  the  system  and  has  the  properties 


The  solution  of  (lL)  satisfying  the  initial  condition 


x (o)  = X 
— ^ o 


IS 


- * + J 5 Ct.x)f  Jr 


(17) 


Let  us  express  $ in  partitioned  form 


u . v 
w Y 


where  each  sub-block  is  a matrix. 

There  are  two  cases  of  special  interest.  First,  suppose  £-0  . 
Then,  the  propagator  functions  (U)  become 

♦ vtfr.Of  j 

Ir  + * Y £ 


The  Jacobian  determinant  (13)  becomes 


J * At  , 


and  the  solution  of  the  transcendental  equation  (12),  which  is  now 


linear,  becomes 


£■*  = v ] 


Thus,  we  reduce  the  formal  solution  to  the  form 


Gr[wio,t)|  +-Y(o.t ) 2 ] t-  Vio,t)  g jixp  (-7}  (19) 


where 


S j T,  A(V)  It 
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From  this  it  follows  that  the  spatial  density  reduces  to 


! &[«Mk-Y(o.t)V'(o,tJo  (2i) 


where 


M(o,t)  = w <o,t ) - Vio.ti V *’( o,t>  U(o,t) 


(22) 


The  fundamental  matrices  appropriate  for  linearization  about  a 
circular,  Kepler ian  orbit  of  mean  motion  w.  are  (Heard,  1976) 


U(t,T)  r 


W(t/0  = 


/ c 

o 

h 

2(l-t) 

°\ 

-is 

1 

0 

VU,r)  - i 

-Ul-c) 

Hs-llt-c) 

o 

A 

i 0 

1 

l o 

0 

J 

O 

s/ 

1 ->V5 

o 

° \ 

1 6 

0 

w> 

rl 

1 

o 

o 

1 

o 

VU.T)  » 

O 

1 o 

\ ° 

o 

u 

o J 

and 


I 3nCt-T  ) ♦ A-  S 
-ZO-c) 

\ o 


2(l-t) 

s 

o 


0 \ 

o 

D/s  / 


(23) 


(24) 


where 


s * sU\.  t «an,nit-r) 


and 


p = c.) 


(25) 


In  some  applications  of  this  case  it  is  more  appropriate  to 
linearize  about  an  elliptical,  Keplerian  orbit.  The  fundamental 


matrices  for  this  case  can  be  obtained  from  Tschauner  and  Hempel 

(1965). 

The  second  special  case  worth  elaborating  upon  is  the  case  for 
which  A and  t are  constants.  Omitting  the  details,  one  obtains  for 
this  case 


f(i  .*> 


Ui<  \f) 


O 


VO  J 


l Ob) 


where 


a'  f 


07  ) 


107 


and  all  matrices  are  evaluated  at  the  arguments  ( o ,t  ).  This 
special  case  would  be  useful  if  one  desired  to  include  perturbations 
or  thrusting  through  the  vector  t . 

DIRECT  FR0HIJ5M  ^Disintegration  over  a Finite  Interval  of  Time) 

In  this  section  we  generalize  the  direct  problem  treated  in  the 
previous  section  to  allow  the  particles  to  disperse  from  a point 
source  over  an  interval  of  time.  This  generalization  would  be 
Important  when  considering  the  dispensing  of  orbiting  dipoles  as 
in  the  case  of  the  Westford  Needles  ^Morrow,  et  al. , 1061)  and 
when  mode 11 lug  dusty  comet  tails  (Pinson  and  Probstein,  loot?). 

Tlie  phase-space  distribution  function  now  satisfies  the 
differential  equation 

10 


* 


or 


JLf  = S 


(28) 


M + f A = S 
D t 


where  S(<^,t)  represents  a source  of  5 particles  per  unit  time 
created  in  the  volume  element  d^d^  at  and  at  time  t . 

Let  denote  the  solution  of 


£ 4>  - o 


(29) 


which  satisfies  the  initial  condition 


= . (30) 

The  function  <$>(  <^,  £,t|0  is  'the  phase-space  distribution  for 
particles  emitted  at  location  <^(t)  at  time  t . It  is  defined  for 
all  ^ j and  for  t > r 
Now  let 

t 

dr  . (31) 

Then  -f  satisfies  the  differential  equation 
oCf  = Cr<f)  - <^s(t>] 
and  the  initial  condition 

* O ^«r  ^ ( t , ) . 


11 


Thus,  f is  the  phase-space  distribution  function  for  dispersion  from 
a point  source  whose  trajectory  is  4<t)  and  whose  distribution  of 
momenta  are  described  by  wi;) . 

The  goal  is  to  find  an  expression  for 


To  this  end  we  use  the  propagator  functions  Q 
to  write 


Pr<W'lV*>« 


i.p.t  /t.o-  ^sit)J  \ (32) 


where 


J & [|(0,  f<n,t]  J+  , (33) 

t 

and  the  integral  is  to  be  evaluated  along  that  trajectory  which 
passes  through  J $ at  time  t . 

We  are  now  in  a position  to  obtain  j> < -i, p, f>.  Let  us  denote 

j • 

According  to  the  results  of  the  previous  section  we  may  write 

J ui 

where  f./'l.iW)  is  the  solution  of  the  transcendental  equation 


(3U) 


12 


and  J is  the  Jacobian  determinant 


■y  3 ( Q ) * " * ) Q.  ^ 

J = ( 

p,  - 

Finally,  upon  interchanging  the  order  of  integration,  is 

reduced  to  the  quadrature 


= J p Jr  . 


For  a linear  dynamical  system,  the  results  of  this  and  the  pre- 
vious section  may  be  combined  to  reduce  equation  (32)  to 

V*.tU  )-  g[w(U><j.  + Y<r.*)p]  &[uu,t)^  + irj]  *7  {-  T } 


where 


T=  /TrA^)J 


From  this  we  obtain 


r,lttlr)=  |i^V(r.t)|  G [NN(T,t)^-Y(T,t)V'(T,0  <^,U>]  tvr{-r  \ J 


where 


= W-irt)-Y|r,t)V'|Tt)U(r 


To  illustrate  the  application  of  these  results,  let  us  examine  the 
case  of  a source  moving  on  a circular  orbit  in  a gravitational  field 
whose  epicyclic  frequency  is  exactly  twice  that  of  the  orbital 
frequency.  Such  a condition  obtains,  for  example,  in  the  inner  regions 


13 


— " ■ • 


— 


of  galaxies.  Given  a suitably  contrived  momentum  distribution  func- 
tion, a closed  form  solution  can  be  obtained  in  this  case. 

Let  us  restrict  our  attention  to  the  behaviour  of  the  ensemble  in 
the  orbit  plane  of  the  source.  Let  us  further  choose  units  3uch  that 
the  radius  of  the  orbit  of  the  source  and  the  orbital  frequency  are 
unity.  Then  the  motion  of  the  individual  particles  is  determined 
from  the  Hamiltonian 

>€  = i ( p,1,  ♦ f ‘ - i 5 pt  - 2 f 1 

where  the  radial  distance  of  a particle  is  v®  l*{  , the  angular 
coordinate  is  9 ^ , and  the  conjugate  mementa  are  i and 

K • 1 ♦ * C 

From  the  general  expressions  given  in  Heard  (1976)  we  obtain 

_ I f coi  ( t - e ) iM.  i*  - r ) \ 

M(t.r)a  


and 


J -a  i — v1  ( t - r ■)  , 1 ■=.  o 


Therefore,  the  partial  density  becomes 


5 , i(5s^£) 


(38) 


«<>  lh  5 a . U-T)  , t*«ar(t-T>. 


A choice  of  momentum  distribution  function  which  allows  a closed 
form  solution  is 


Ik 


f.  ) = 


1 


(39) 


1 * **• 

No  particular  physical  significance  is  attached  to  this  form  of 
other  than  it  does  emphasize  the  smaller  ejection  velocities  as 
intuition  would  demand.  Its  main  virtue  is  one  of  mathematical  con- 
venience. 

After  substituting  equations  (38)  and  (39)  into  equation  (36) 
we  obtain 


) = 


i 

X^F 77 


drc  ttu/s* 


% 


t 


(1*0) 


where  ^ and  where  we  have  taken  t0  ■=  o . In  using  equation 

(1*0),  care  must  be  taken  to  choose  the  proper  branch  of  the  arctangent 
function. 

A particularly  revealing  form  for  may  be  obtained  for 

large  values  of  ^ , say  <^  >>  1 . In  this  case,  an  expansion  of 
may  be  obtained  in  powers  of  1 J y- 
We  have 


J 


Hr 


1 ‘ ♦ ■ **>■*  \ t - r ) 


t 


Therefore 


fl\  A Vj*  - (V2]0(t-  i * ° ( vy  > 


(>*i) 
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This  shows  that  the  density  is  asymptotically  linear  in  time  with  a 
small  periodic  fluctuation. 

A serious  physical  objection  to  the  choice  of  equation  (39)  for 
the  momentum  distribution  function  is  that  there  is  no  upper  bound  to 
the  velocities  of  the  ejecta.  Consequently,  as  soon  as  ejection 
begins,  particles  appear  at  all  locations  in  space.  To  remedy  this 
defect,  we  may  modify  equation  (39)  to  require  an  upper  limit  cut-off 
for  the  velocities  as 

C(e'  * TT^-  H < r.1  - c">  • (>*2) 

The  resulting  expression  for  o^.t)  is: 


8 


where  S is  the  union  of  intervals  over  which  the  inequality 

P„‘  > <\  / (43) 

is  satisfied. 

For  o «.  t c rr  we  have 

8 = 0 

-[S.t] 


, t S = Arc  ri.  ( ^ / f>.) 
, t > 6 


1 6 


For  tt  t ^ Ltt  we  have 


C>  = £ S , IT  - & ] 

- [ S,  Tf-  s J U [n-t-Sjt] 


, t •<-  ir  *•  8 

> t > it  ♦ 8 


and  so  on  . 

An  example  of  numerical  evaluation  of  the  above  is  shown  in 
Figures  (1-2).  A descriptive  account  of  the  evolution  is  as  follows. 
The  ensemble  is  always  symmetric  about  the  origin,  i.e.,  bounded  by 
a circle.  It  is  originally  concentrated  at  the  origin.  It  grows  in 
size  until  it  reaches  its  maximum  radius  of  p.  at  i = v . Thereafter 
the  density  increases  everywhere  in  the  circle  as  particles  are 
continuously  ejected  but  the  radius  of  the  circle  enclosing  the 
ensemble  remains  p0  . 

An  example  of  wider  applicability  is  the  case  of  continuous 
dispersion  from  a circular,  Keplerian  orbit.  The  quadrature  (37)  must 
be  performed  numerically  in  this  case.  An  example  of  such  an  evalua- 
tion is  shown  in  Figure  3.  These  results  are  based  on  an  isotropic, 
Gaussian  momentum  distribution  function  and  show  isodensity  contours 
when  emission  has  occurred  over  one-quarter  of  a revolution  of  the 
parent  body. 

THE  ASYMPTOTIC  PROBLEM 

The  purpose  of  this  section  is  to  describe  the  asymptotic  state 
( in  the  limit  t — m ) of  the  ensemble  when  the  particles  emanate  from 
a point  and  when  the  Hamiltonian  for  the  individual  particles  is 


/,  I 
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Fig.  1 - Spatial  density  versus  q at  various  times 


contours  at  1 3 »/2  for  a continuous  source  in  a 
Circular,  Keplenan  Orbit  (Spherical  Velocity  Distribution) 


20 


separable.  The  asymptotic  behavior  of  the  distribution  function  is 


determined  by  applying  theorems  developed  by  R.  T.  Prosser  (1969)  in 
his  spectral  analysis  of  central  force  motion.  The  theory  leads  to 
the  weak  limit  of  the  spatial  density  function  and,  implicitly,  to  the 
boundary  of  the  asymptotic  domain.  The  theory  is  illustrated  by 
application  to  the  cases  of  elliptical  Keplerian  motion  and  to  the 
case  of  motion  about  an  oblate  primary. 

Kotsakis  (ldoU)  has  studied  the  case  of  Keplerian  motion  for 
isotropic  dispersion  from  a circular  parent  orbit.  He  obtains  the 
boundary  of  the  domain  as  the  envelope  of  a family  of  ellipses.  He 
also  obtains  a radially  averaged  density  of  the  particles.  The 
results  from  our  dynamical  theory  generalize  this  to  anisotropic 
dispersion  from  an  elliptical  orbit  and  provide  the  radial  dependence 
of  the  spatial  density.  Wiesel  (1976)  considers  the  related  problem 
of  a statistical  theory  of  the  Kirkwood  gaps.  The  basic  difference 
between  Wiesel' s theory  and  the  present  one  is  that  we  invoke  Prosser's 
theorems  where  he  uses  ergodicity  on  integral  surfaces.  Whereas  our 
results  are  strictly  limited  to  separable  systems,  Wiesel '3  theory  has 
the  potential  of  wider  applicability. 

Mathematical  Formulation 

Consider  a dynamical  system  defined  by  a separable  Hamiltonian,  X.  . 
The  trajectories  of  this  dynamical  system  define  a flow,  T,  ; on  phase- 
space  by 

rt  : j < r ) *■  v , t - t a 

i „ 

where  x \ » i . ; and  ^ . » are  the  coordinates  and  conjugate  momenta. 


I 


where  i > \ denotes  the  Poisson  bracket.  Let  f u>  = F(*.t> 

and  ^ - “(x,o)  . The  formal  solution  of  the  partial  differential 

equation  (UL)  may  be  written  in  the  above  notation  as 


Let  S be  a canonical  transformation  from  * to  some  linear  com- 
bination j.  , w of  action  and  angle  variables, 

S ; j..* 

and  let 


be  the  fundamental  frequencies  of  the  motion.  The  composition  of  the 
flow  mapping  and  the  canonical  transformation  yields  the  mapping 

STt  : < ^ f ) > ( j.  , 'i  v t v ) . 


For  particles  which  disperse  from  a common  point,  the  initial 
phase-space  distribution  function  may  be  written 


The  first  factor  is  the  function  describing  the  initial  distribution 


of  the  momenta  and  is  assumed  to  be  known.  The  second  factor  is  the 


Dirac  delta  function  and  it  describes  the  fact  that  the  particles 


originate  at  coordinates  . 

Since  the  transformation  S is  canonical,  and  therefore  has 


Jacobian  determinant  unity,  we 


where  the  transformation  S is  expressed  in  the  form 


*The  application  of  the  theorems  of  the  previous  section  to  generalized 
functions  has  not  been  Justified  rigorously  here.  To  do  so,  one  would 
regard  a generalized  function  as  an  equivalence  class  of  regular 
sequences  of  good  functions  (Jones,  1966)  and  develop  arguments  for 
the  validity  of  the  interchange  of  limits  and  integrals  . 


Define  the  vector  valued  function  w,  implicitly  by 


9 1 j .*.■)  = 


so  that 


Then,  from  a fundamental  theorem  about  the  change  of  variable  for  Dirac 
delta  function  (Jones,  19bb)  we  may  rewrite  (!*7)  as 


where 


£ U ] j 


9 w(>,w) 


If  equation  (U8)  has  more  than  one  root,  the  right-hand-side  of 
equation  (1*9)  must  be  evaluated  at  each  root  and  summed  over  all  of 
them. 

We  must  allow  for  the  possibility  of  degenerate  systems,  such 
as  Keplerian  motion,  in  which  some  of  the  frequencies  are  identically 
equal  to  zero.  To  do  this,  let  us  partition  the  variables  into  two 
parts 


2S 


J 


\ w*"/ 


» * 


} 


such  that  the  frequencies  3/-  X , auid  satisfy 

t ^ 

K =0  l«..t  ) 


and 


0 


(*■■*■) 


Then,  from  Corollary  1,  the  weak  limit  of  the  phase-space 
distribution  function  is  expressed  in  new  variables  as 

f 4,w  ) - S ['«£'- w:«h]/j,(  f)  • 

In  terms  of  the  original  variables  the  result  is 

f{s[n,.»)]}  Hr., )A[j.V(.j  <5°> 

The  asymptotic  form  of  the  spatial  density  function,  Jin  , 
is  obtained  by  integrating  the  f of  equation  (50)  over  all  momenta 
S . There  are  two  cases  to  be  considered.  First,  if  the  set  w 
is  vacuous  (i.e.,  the  system  is  non-degenerate)  we  obtain 


(51) 


Second,  if  the  set  w'  is  non-vacuous  (i.e.,  the  system  is  degenerate) 
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(52) 


FCV>-  Jjj(vid)  ^ 

E 

where  the  integration  is  performed  over  the  hypersurface 
E » W'c^.p)  -*.[  J(^,*)]  - O , 

and 


r = o 


Two-Dimensional  Keplerian  Motion 

In  this  section  the  details  of  the  algorithm  are  worked  out  for 
the  case  in  which  the  particles  execute  bounded,  Keplerian  motion  in 
a common  plane.  There  are  three  reasons  for  considering  this  example 
in  detail.  First,  it  is  simple  enough  that  the  algebraic  details 
do  not  obscure  the  structure  of  the  algorithm.  Second,  from  the  work 
of  Kotsakis  (196H)  we  have  a partial  check  of  the  algorithm  for  the 
more  specialized  case  of  isotropic  dispersion  from  a circular  orbit. 
Third,  when  the  velocity  increments  are  small  relative  to  the  total 
velocity  in  the  parent  orbit,  the  out-of-plane  component  is  de- 
coupled from  the  in-plane  components  (Kotsakis  196U,  Tschauner  and 
Hempel,  1965).  Thus,  this  case  is  really  a prerequisite  for  the  full, 
three-dimensional  Keplerian  case. 

The  analysis  will  be  executed  in  three  stages.  First,  the 
expression  for  J for  a general  momentum  distribution  function  P 
will  be  reduced  to  a one  dimensional  quadrature  which  normally  will 
have  to  be  performed  numerically.  Second,  the  quadrature  will  be 
performed  analytically  for  the  case  of  anisotropic,  small  velocity 
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dispersion  from  an  elliptical  orbit  when  ? may  be  represented  as  a 
Dirac  delta  function.  Third,  the  results  will  be  compared  to  those  of 
Kotsakis  and  numerical  examples  will  be  presented. 

A.  Reduction  to  Quadrature 


Let  the  generalized  coordinates  be  the  polar  coordinates 


whose  conjugate  mementa  are  £ •=  \ r‘  e / • 


The  Hamiltonian  is 


= i ( K * K ) ~ ^ /*■ 


Let  the  coordinate  system  be  oriented  such  that  the  initial  point  is 


The  Delaunay  variables  are  a convenient  combination  of  action  and 


angle  variables  so  w and  j.  become 


«•(!) 


» V 3 


The  system  is  degenerate  so  the  angle  variables  and  frequencies  are 


partitioned  according  to 


v'=  o 

v " - /UL  ‘ L 


The  well-known  relations  connecting  ( ^.  p ) and  ( w i ) are 


L = /*■  ( 1/>/r  - P,1’  - pvVr  «•  ) 
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where  e'  _ I - g/l1-  and  v(e,2)  is  a periodic  function  of  2.  . 
The  initial  phase-space  distribution  function  is 

f.  ~ S(  r-r.1)  6(0  i/ra  . 


r 


The  remaining  task  is  to  obtain  J which,  as  we  recall,  is 


J - 1 $ J P.  J N • 


(55) 


To  evaluate  the  integral  (55)  it  is  convenient  to  transform  coordinates 
from  p,  , ^ to  n ^ . In  so  doing,  we  obtain  the  one-dimensional 
integral 


/ = J J !V 


111 

a» 


.El. 


IJI 


(56) 


The  expression  for  the  argument  of  $ is 


nfp,<o,r<),Pt]a  -~^(i  i().,fc£))l  (57) 

s*«.  <3  A*-  r / p 


and 


Equation  (57)  is  most  directly  derived  by  eliminating  ^ from 

r“  P*  />*•[  ' + * 
r°  ^ ^ ( I *•  * ““  3 ) 


then  solving  for  pi  as  a function  of  , and  finally  using  the  energy 
integral  in  the  form 

*•  NV*-1-  - ly*/r  - \1  + |\/r/  - =>A,  * 

This  derivation  allows  us  to  interpret  ^ ft')  as  the  radial  velocity 
which  sends  an  orbit  with  initial  conditions  *v  3 j t , ft.  through  the 
point  r,  d . 

i 

■*«aa 


The  third  factor  in  the  integrand  of  (56)  involves  ip,/?*]  which 


!Ll f (21®  .l)  + 

?*i  s^9[r  \ r»  r ) pa 


Finally,  we  assemble  the  above  results  to  obtain  the  desired  one- 
dimensional integral  for  o : 


f-m 


'T  OF 


where 

fi  ~ P» />»*•”  1 * /30  = f\//L*-’r.  - 1 , c:w9  * - s*~  9 • 

This  solution  is  valid  for  general  momentum  distribution  functions 
■^lPt,Pi')  • However,  in  most  cases  it  would  be  necessary  to  resort  to 
numerical  methods  to  evaluate  the  integral.  In  the  next  section  we 
consider  an  important  case  for  which  the  integrand  (59)  can  be  evaluated 
analytically. 

Small  Initial  Velocity  Increments  Concentrated  on  an  Ellipse 


Let  us  now  consider  the  case  where  the  momentum  distribution 


function  is 


This  describes  anisotropic  dispersion  from  the  parent  orbit  whose 
initial  conditions  are  r, ? 0 •’  f'.g  , Po  • The  anisotropy  is  of  the  special 
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form  that  the  incremental  momenta  of  the  particles  lie  on  an  ellipse 
in  the  p, , -plane.  Were  <r,  = <rz  , the  dispersion  would  he 

isotropic  and  the  incremental  velocity  of  each  particle  would  he  j-,  v0  . 
Were  P,,=  o and  p*/^  = r.  » the  parent  orhit  would  he  circular. 

The  integral  (59)  may  now  he  written 


/=  J I «■***.* 


3/1 


™ Vf  (Tt  fx 


(6l) 


where 


f ( 


M = 


g0c-/3 


~ Pm 


+ /AIM1  - v* 


O1 


The  integral  (6l)  may  he  evaluated  immediately  as 


}(*.*)*  y JL_i 

where  p*  is  the  root  of  the  equation 

? < p * ) = O . 


3/1 


(62) 


p**p* 


(63) 


Equation  (63)  represents  a quartic  which  has  two  real  roots,  and  the 
summation  in  equation  (62)  is  taken  over  both  of  these  roots. 

In  most  cases  of  physical  interest  r.  v8/p0«1  , meaning  that  the 
incremental  velocities  are  much  smaller  than  the  total  velocity  of  the 
parent  body.  We  shall  exploit  this  fact  to  solve  the  quartic  (63) 
and  to  give  the  results  in  their  most  illuminating  form. 

To  this  end,  we  let 
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a)  = 


T.Vm  c - l 


F.< 


^ ] 


(61*) 


r s i p«  H ” ♦ * ) 


and 


P*  = K i 1 ♦ « (*)  ] 


where  !*|  , |t|  «.  7 — 1 


Mo)  = o 


(65) 

(66) 


The  curve  rt9)  9)  is  the  orbit  of  the  parent.  Upon  substitu- 

tion of  equations  (65)  - (66)  into  equations  (62)  - (63),  and  after 
retention  of  only  the  lowest  order  terms  in  v.  , one  obtains 


J ( * , 9 ) ~ 


n a‘.'T»  Fj  N.  m, 


A 


where 


(67) 


= 


,=(ir) 


t. 


(68) 


In  equation  (67)  we  have  used  the  relations 
- P,0  P*/>  = e.  j. 

where  e.  and  are  elements  of  the  parent  orbit. 

.he  curve  bounding  the  asymptotic  domain  is  obtained  by  setting 


oo 


I I 


O . Thus,  the  bounding  curve  is  defined  by  the  equation 


X.  „ 19} 


p.  I '/*»■»'  ‘ 


(69) 


C.  Discussion  of  the  Two-Dimensional  Case 

First  we  should  take  the  opportunity  to  compare  our  results 
to  those  obtained  by  Ketsakis  in  the  special  case  of  isotropic 
dispersion  from  a circular  parent  orbit.  If  we  let  p .‘,^^  43,  and 
r,  rj;=i  ; equations  (69)  and  (63)  become,  respectively, 


=■  -a  V<./r  v0  [ ] 


(70) 


or 


and 


or 


x.  t i-3(l  + 3 s**-1  1 9 "> 

f » ■ - X1  j 

J Tf  D * ou  W 


Vl 


_ »/2 


Y •=.  — - — - { * ^ \ 1 * * S¥tK*  ^ ] 


(71) 


where  «.  = v./v,  and  vt  - /*»/*•. 

Equation  (70)  is  equivalent  to  equation  (lU)  of  Kotsakis. 
However,  Kotsakis  does  not  derive  J^.d)  but  rather  the  radially 
averaged  density 

♦ 

J ? 

p Id)  = 


r - U 

2+X^  J J j3v*ir^9 

* x. 


3k 


i 


***•  -J* 


I 


Using  equation  (71)  we  calculate 


■=.  l/ ( i-~  r/  ) (4<x  ) j — . 1 a ( 1 ♦ 3 1 9 ) 


~ ’ / *♦  r.T,  fl~-  t 0 ( 1 ■*  j 1 9 ) ' 1 


which  agrees  with  equation  (21')  of  Kotsakis. 

The  shape  of  the  bounding  curve  is  quite  sensitive  to  the  value 
ol  when  e „ is  not  small  and  there  is  also  a marked  dependence 
on  the  relative  values  of  <r,  and  crt  . Figures  (Ua-c)  show  these 
bounding  curves  for  the  case  e0  = o.v  and  for  various  values  of  ^ 
and  . 

The  radial  dependence  of  the  density  is  illustrated  in  Figure  (5). 
The  salient  feature  is  the  singularity  at  the  inner  and  outer  ex- 
tremeties  of  the  domain.  These  singularities  are  artifacts  of  the 
singular  nature  of  the  function  chosen  for  £ and  would  be  smoothed 

out  upon  a superposition  of  a continuum  of  such  functions. 

It  is  interesting  to  invert  the  problem.  That  is,  to  suppose 
that  one  observes  an  asymptotic  distribution  of  particles,  known  a 
priori  to  result  from  a momentum  distribution  of  the  form  (60).  and  to 
ask  if  we  can  determine  this  momentum  distribution  function.  Clearly, 

this  involves  only  the  determination  of  the  parameters  p,0)p.  r,  9,/,,^, 

and  v,  . If  one  knows  the  boundary  of  the  domain  one  can  determine 

r O),  ra  and  9,  and  hence  Pw  and  simply  by  tracing  the 

center  of  the  domain.  Further,  the  shape  of  the  boundary  yields  v, 

and  v»<r.  unambiguously.  To  separate  v,  from  these  latter  paramters. 

however,  it  is  necessary  to  know,  in  addition,  the  density  at  an 
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Fig.  4a  - Boundary  for  the  Keplerian  asymptotic  domain 
in  the  orbit  plane  of  the  parent  body  for  e m 0.7  and 
various  values  of ~r.  and  -r-.  The  dotted  curve  is  the 
orbit  of  the  parent  body.  Here  g « 0°. 
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interior  point.  The  central  density  at  any  point  for  which  <9  * 
would  suffice.  We  conclude  that  knowledge  of  the  boundary  and  the 
density  at  one  interior  point  are  sufficient  to  solve  the  inverse 
problem.  Such  complete  knowledge  is  not  necessary,  however.  For  if 
we  know  the  point  of  origin  and  the  boundaries  at  two  more  distinct 
values  of  9 we  can  unambiguously  fix  r. , , Pl0  > p„  .''.a-, , and  v„ 

Knowledge  of  the  density  at  an  interior  point  would  then  be  sufficient 
to  solve  the  inverse  problem.  On  the  other  hand,  even  complete  know- 
ledge of  the  density  along  just  one  ray  9 a const,  is  insufficient 
to  solve  the  inverse  problem. 

Three-Dimensional  Keplerian  Motion 

In  this  section  we  provide  the  details  of  the  algorithm  for  the 
case  of  general,  bounded  Keplerian  motion.  As  in  the  previous  section, 
the  general  case  will  be  reduced  to  a one-dimensional  integral  and  then 
the  quadrature  will  be  executed  analytically  for  the  case  of  small 
velocity  disperison  and  incremental  velocities  confined  to  an  ellipsoid. 

A.  Reduction  to  Quadrature 

Now  let  the  generalized  coordinates  be  spherical  coordinates 

whose  conjugate  momenta  are 

* = 

The  Hamiltonian  is 

X = J < P,L  * PvV/r  ‘ * 'l  - AA./r 


OI  ‘ 0 0 


!*0 


“ — ’c-"'  n"  ' 


The  initial  phase-space  distribution  function  is 


~ ( P\  . ft , pj  ) S ( r - r.  ) S i9)  S ».  f V r,*1 


Using  the  relations  (73)  and  equation  (50),  one  obtains  the  follow- 
ing limit  of  the  phase-space  distribution  function 


y r,,  Pt.ft)  = 5(i)/i 


where 


*.  < F* , P,  ) = <•*«•  ft  . 


^ V(4  -l)j  -a«,.  j^-  V 


£ * •P  ~ a,‘  J»~  ^ ( 0 j 


In  the  last  two  equations,  the  variables  L,  G,  H must  be  regarded  as 
functions  of  J>t  f,  r and  9 as  given  in  equations  (73).  In 
lerivtng  these  expressions,  use  has  been  made  of  the  fact  that  ^ 
and  the  initial  value  of  ».  are  either  O or  it  because  of  the  choice 
of  the  o-Tdinate  system.  The  Jacobian  determinant  appearing  in 


equation  CM  is 


J.  - Hi.  j.  k) 


1*2 


The  limiting  spatial  density  function  is  now 

Fs  If  • ( } 

To  evaluate  the  integral  (75)  it  is  convenient  to  transform  coordinates 
from  to  7,^  f>,  to  obtain 


f = J S'WftfoA , fi)<  pj  } 

X 9 C f » ■ PW  ) I <tp; 


Omitting  the  details  of  the  derivation,  the  final  expressions  for  the 
first  two  arguments  of  S in  equation  (76)  are 

• • * • 

*1 1 M = r,).ftfo.0.i>3 ) t ^ ] 

= ^ M <A-£»  4 tar  9 ; *•  ( YY  J 

pj  ( 1 + Si*\x9  1 /6  J I — <«>.*  • 9 *-o:  1 j6 

and  = 'C  [ p, ) j |»j  3 

= P;  Wh's/.'^/s  ^ (l8) 

where  - ( pf/ru. ) ( ® v 9 / fw,1  ■}>  ) — 1 


f 


"p,",w"-,-pl '*■  ,,J|  1 '•  


and 


4o  * j.‘«c‘9  * +<w*  9 / ^ J - 1 

For  the  second  Jacobian  determinant  in  equation  (76),  one  obtains 


f.,K>  _ 


3(  1 , t ) 


(79) 


( 1+  io~v9  O-  ) 

where  it  is  understood  that  the  Jacobian  has  been  evaluated  for 

Assembling  the  above  results,  we  obtain  the  following  one-dimen- 
sional integral  for  f : 


f = J i 'A-A.™*  h 


/ab  fj 


, 3/v 


/ 1 S-yfly  (80) 

-os9  i'*~?  J (AS)3'1"  r* 


where  the  abbreviations 

A = ( 1 + :1~*"  © to+  l^>  ~)  / co: 1 3 

5 S 1 — CoS  1 ^ Co."  * $ 

have  been  introduced. 

Equation  (80)  reduces  the  calculation  of  J>)  for  arbitrary 

> to  a single  quadrature.  As  in  the  preceeding  section  we 
next  consider  an  important  special  case  for  which  the  integral  (80) 
can  be  evaluated  analytically. 

B.  Small  Initial  Velocity  Increments  Concentrated  on  an  Ellipsoid 
Let  us  now  consider  the  case  where  the  momentum  distribution 


function  is 


■■ 


V 


- 'f, 


1 


( 8l) 


i 


!> 


r.'-o-x 


<■  Pi  - p*  >' 


This  describes  anisotropic  dispersion  from  a parent  orbit  whose  initial 
conditions  are  r*,0  . The  anisotropy  is  of  the  special  form 

that  the  momentum  increments  lie  on  a triaxial  ellipsoid.  Were 

= , the  dispersion  would  be  isotropic  and  the  incremental 

velocity  of  each  particle  would  be  <n  v„  . Were  plQ=  J and 
the  parent  orbit  would  be  circular. 

The  integral  (59)  becomes 


where 


H_  ii.’-l  6 -mi  -P  + -pa.  co.9  p | 

{cp  )l  " > (82) 

r J p*  > ^Ai3  )'  ( rr<r,  r*  V,1  ) 


f\>  = 


1 

" 2 

1 i <3  -/3„  co;  9 S yj-g- 

- 

P,  9 

AftO 

P.j  + j 

P,  J 1 

< r,  - k r _ 

I 


Now,  we  again  exploit  the  fact  that  the  incremental  velocities  are 
expected  to  be  small  relative  to  the  velocity  of  the  parent  body, 
and  write 


r - ( f’W/u.  ) ( r v «.  ) 

Pj"  = < i ♦ * ) 

I 

where  r is  given  by  equation  (6U);  i*.i,i«I  and  I9l«-l  ; and  f k 

U5 

kt_- 


i 


’ 0 . 


It  follows  that  A « 1 -t-0( ®l)  and  Bs  ^ , After  retaining  the 

lowest  order  terms  in  x and  ® in  the  evaluation  of  the  integral  (82) 
one  obtains 


a"'/1 


J ^ <rv  <r,  rQ  v„: 


where 


s •?  s 4 t c * 'ot 


A , [(  r Jl  l /r^('  _©J_  . (8U) 

/ *1  <r, 1 J _v  p*  ) - a-, r **■ 


The  asymptotic  domain  is  bounded  by  the  surface  defined  by  A = o , 
that  is,  the  surface 


? )*■ 


9 )i~ 


- v1  (85) 


m 


c.  Discussion  of  the  Three-Dimensional  Case 

Following  the  procedure  of  the  previous  section,  it  is  easily 
shown  that  equation  (8U)  reduces  to  that  obtained  by  Kotsakis  in  the 
case  of  isotropic  dispersion  from  a circular  orbit.  Similarly,  if 


where  a and  b are  the  semi-axes  of  the  cross-sections  of  the  domain  by 


const.,  one  agrees  with  the  averaged  density  reported  by 


Kotsakis  except  for  a factor  of  two.  This  discrepancy  arises  because 


of  an  apparent  slip,  by  a factor  of  two,  in  Kotsakis'  calculation  of 


The  cross-sections  of  the  asymptotic  domain  by  a plane  const 


remain  elliptical  as  in  the  isotropic,  circular  parent  orbit  case 


In  sections  of  the  bounding  volume  by  the  planes  £ const 


the  iso-density  contours  are  ellipses.  The  density  is  again  singular 


at  the  boundary  of  the  domain.  The  density  minimum  is  attained  on  the 


parent  orbit 


The  cross-section  of  the  bounding  domain  by  the  plane  6 ~ -> 


is  the  bounding  curves  described  in  the  previous  section 


one  must  replace  <j\.  by  <r, 


Particles  Orbiting  an  Oblate  Primar: 


The  results  obtained  in  the  previous  section  are  unrealistic  in 


the  sense  that  perturbations  due  to  the  oblateness  of  the  primary 


or  the  attraction  of  a third  body,  would  precess  the  orbits  and  destroy 


the  highly  stylized  structure  of  the  Keplerian  case.  One  might  argue 


that  when  the  perturbations  are  small  there  will  be  an  intermediate 


period  of  time  when  the  Keplerian  dispersion  will  dominate.  During 


this  intermediate  period,  the  structure  we  have  just  discussed  would 


provide  a decent  approximation  to  the  state  of  the  ensemble.  Never 


theless,  in  the  limit  t 


the  perturbations  would  eventually 


dominate 


Fortunately,  the  theory  developed  in  Section  3 is  sufficiently 


we  shall  use  the  algorithm  to  determine  the  boundary  of  the  asymptotic 


domain  when  the  particles  orbit  an  oblate  primary.  The  discussion 


of  the  density  in  this  case  will  be  left  for  a future  paper.  In 


order  to  apply  the  theory,  it  is  necessary  to  deal  with  systems 


governed  by  a separable  Hamiltonian.  The  Hamiltonian  describing 


motion  about  an  oblate  primary  is  not  separable,  except  in  the  case 


of  the  Vinti  potential.  However,  the  motion  of  a particle  about 


an  oblate  primary  may  be  described  accurately  by  an  intermediary  orbit 
which  is  based  on  a separable  Hamiltonian  (Garfinkel,  1959,  Aksnes, 
1972). 

A.  Boundary  of  the  Asymptotic  Domain 

Let  us  again  choose  spherical  coordinates  to  describe  the 
motion  of  the  particles.  Let  the  plane  9 - 0 be  the  equatorial  plane 
of  the  oblate  primary.  Let  the  point  , 0)  be  the  initial  point. 

The  Hamiltonian  for  a spherical  coordinate  intermediary  orbit 
is  (Aksnes,  1972) 


( P-I1"  «"  ri’/r’-'t-  yf/r'-iaS-  e ) - A*-/r  Px  < ) j (86) 

where  is  the  second  zonal  geopotential  coefficient  (~10  for  the 
earth),  is  the  equatorial  radius  of  the  primary,  is  the  second 
Legendre  polynomial,  and  e,  is  a constant  which  may  depend  on  the 
momenta.  Associated  with  this  Hamiltonian  are  Delaunay  type  variables 
which  we  shall  use  as  a convenient  combination  of  action  and  angle 
variables.  The  relations  connecting  these  generalized  Delaunay 
variables  to  the  coordinates  and  momenta  are  more  complicated  than 
equations  (73).  For  our  purposes  we  need  list  only  the  relations 

L = /*■  - p,1-  - )‘,,V  ) 

G-  - ( Pi  ¥ / IOJ*0  ) & - ) ■=.  0 t 

H - (87) 


L 


The  importance  of  the  Hamiltonian  (86),  in  the  present  application,  is 
that  it  is  separable,  that  the  frequencies  associated  with  L,  G,  H are 
non-zero  and  linearly  independent  almost  everywhere,  and  that  it 
accurately  describes  the  motion  of  a particle  about  an  oblate  primary. 

Since  the  dynamical  system  is  non-degenerate,  it  is  not  necessary 
to  partition  the  angle  variables,  and  the  integral  (51)  gives  the 
spatial  density.  This  integral  becomes 

(88) 

where 

4 &.1-  (Vrl  - \/r*  ) +■  -Mr)  } 

S G • — Yl  / la'  1 C | ^4  ) 

and 

G = • 

We  see  immediately  that  coordinate  is  absent  in  equation  (88  ) so 
that  the  density  and  the  bounding  surface  are  rotationally  symmetric 
about  the  axis  of  symmetry  of  the  primary. 

As  long  as  J ^ is  not  zero,  no  matter  how  small,  we  must  deal 
with  the  three-dimensional  integral  (88)  instead  of  the  one-dimensional 
integral  derived  for  the  degenerate,  Keplerian  case.  Having  accounted 
for  this,  however,  we  may  now  neglect  the  effects  of  J0  and  still 
obtain  an  accurate  description  of  the  state  of  the  ensemble.  A 
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consequence  of  this  is  that  the  relations  (87)  may  he  replaced  by 
the  relations  (73)  in  what  follows. 


To  discuss  the  boundary  of  the  asymptotic  domain,  we  shall 


treat  the  following  form  of  y\p> 


where 


are  the  initial  conditions  for  a circular 


parent  orbit.  The  current  choice  of  coordinates  does  not  allow  us 


to  assume  that  as  in  equation  (8l).  From  equation  (89)  it 


follows  that 


where 


In  equation  (91),  use  has  been  made  of  the  relations 


where  i is  the  inclination  of  the  parent  orbit. 

A point  (r  9 ^ > lies  in  the  asymptotic  domain  of  the  ensemble 
if  the  equation 


5 ( r q 


p,  . P-v  , Pm  ) - o 


has  a solution.  This  equation  alone,  however,  is  insufficient  to 
determine  the  domain.  We  must  also  satisfy  the  constraints 


- f,'1'  ( ) ->  o 


p,*"  ''•Yx  *■  pj~  K l/r1-  - I />,*■  > + - t/r  ) > O . 


Now  let  us  concentrate  on  the  case  of  smad.1  velocity  increments. 


This  allows  us  to  represent 


and  to  assume 


oc  ~ QC  x.  ) t I m.  V 1 


Further,  let  us  choose  units  such  that  /*>  - I <\  =■  I . With  these 
simplifications,  equations  (93)  - (95)  become 

f’*-  ♦ !>«.  - (\»  n ]*■  + c pj  - P,4>  I ♦ J*-  Jx  , 


Jt’--*  !».«• 


* 3 »*- 


<■  «■  t3kl) 
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p,w  fj  iui9.)  (jx1-;*)  > *.(*“-*) 


vp*  > p,v  ( - -««.'■£«  ■) 


Equations  (96)  - (98)  provide  a useful  geometrical  interpreta- 
tion of  the  problem  of  locating  the  boundaries  of  the  asymptotic 
domain.  Regard  equations  (96)  - (98)  in  coordinates  p,  _ rv  p., 

Equation  (96)  represents  an  oblate  spheroid  centered  at  [ o,  ' «■  O 
] . The  boundary  of  the  region  for  which  equation  (97) 

is  satisfied  is  either  a triaxial  ellipsoid  ( x < 0)  or  an  elliptical 
hyperboloid  of  one  sheet  (x  > 0),  both  centered  at  the  origin.  If 
9 < , equation  (98)  is  always  satisfied.  If  , the  boundary 

of  the  region  for  which  equation  (98)  is  satisfied  is  a wedge.  Thus, 
our  problem  is  to  find  those  pairs  (x,  9 ) for  which  the  ellipsoid 
(96)  intersects  or  lies  within  the  wedge  (98)  and  the  ellipsoid/ 
hyperboloid  (97). 


£ m S oc  y/coi  4 - .*>*  X t / » 


When  l®l*i  the  equation  and  the  constraints  are  satisfied  if 


<.  1 w «- 


(100) 


This  part  of  the  asymptotic  domain  is  independent  of  9 and  is  the 
volume  between  the  spheres  riltH  , but  outside  the  cones  9 = t.  I 

When  , the  equation  and  the  constraints  are  satisfied 

if 


where 


and 


0 = i v t 


= ><•  [ 1 - (•'  ^ ] 


(101) 


(102) 


The  asymptotic  domain,  as  described  by  equations  (100)  and  (101), 
is  illustrated  in  Figure  7.  This  figure  shows  the  cross-section  of 
the  domain  by  a plane  4>  * const.  The  domain  is  the  volume  of 
revolution  swept  out  as  the  curve  is  rotated  about  the  axis  of 
symmetry  of  the  primary.  An  important  feature  of  the  bounding  curve 
is  that  it  depends  on  the  initial  latitude.  The  perturbations  do 
not  destroy  the  record  of  the  initial  point.  In  the  FCeplerian  case, 
the  asymptotic  domain  is  bounded  by  a toroid  which  is  aligned  with 
the  plane  of  the  parent  orbit  but  is  constricted  to  a point  at  the 
place  of  origin  and  pinched  to  a line  180°  away  in  true  anomaly. 

When  oblateness  perturbations  are  considered,  the  asymptotic  domain 
is  a toroid  whose  axis  is  the  axis  of  symmetry  of  the  oblate  primary. 
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THE  DIVERSE  PROBLEM 

There  is  a large  class  of  problems  in  the  physical  sciences, 
generally  known  as  inverse  problems , which  are  concerned  with  the 
determination  of  models  of  physical  systems  from  indirect , and  usually 
incomplete,  observational  data.  These  problems  abound  in  geophysics 
(Parker,  1977;  Gilbert,  1971)  and  one  may  cite  as  examples  the 
problems  of  determining  the  density  of  the  earth's  interior  from 
seismic  data,  the  structure  of  the  ionosphere  from  radio  wave  propaga- 
tion characteristics,  and  the  structure  of  the  geomagnetic  field  from 
satellite-born  magnetometer  measurements.  Excellent  expositions  of 
the  mathematical  aspects  of  these  problems  have  been  given  by  Keller 
(1976),  Prosser  (1977),  and  Backus  (1971). 

In  a sense,  the  problem  of  orbit  determination  is  an  inverse 
problem.  Certainly  there  sure  similarities  between  the  techniques 
used  by  Backus  and  Gilbert  to  solve  inverse  problems  and  the  differential 
correction  schemes  used  to  estimate  orbital  elements  from  observations 
of  the  position  of  a celestial  body.  However,  there  is  also  a funda- 
mental difference  between  estimation  problems  and  inverse  problems . 

The  former  involve  the  estimation  of  a finite  number  of  parameters 
from  an  over  determined  data  set  while  the  latter  seeks  a continuous 
model  from  a woefully  inadequate  data  set.  It  is  useful  to  draw  the 
analogy  between  the  relation  of  maximum-minima  problems  of  elementary 
calculus  to  problems  of  the  calculus  of  variations,  one  one  hand,  and 
estimation  problems  to  inverse  problems  on  the  other.  A good  example 
of  an  inverse  problem  in  modern  day  celestial  mechanics  is  the 


determination  of  atmospheric  density  models  from  the  study  of  decaying 
artificial  earth  satellites. 

In  this  section  we  consider  the  inverse  problem  associated  with 
an  ensemble  of  particles  which  emanate  simultaneously  from  a point 
source  with  a continuous  distribution  of  momenta.  The  problem  is  to 
determine  the  initial  momentum  distribution  from  observations  of  the 
evolving  ensemble. 

The  inverse  problem  for  the  particle  ensemble  system,  or  any 
system  for  that  matter,  is  completely  posed  only  when  one  has  specified 
precisely  what  observations  are  to  be  used.  For  an  ensemble  of 
particles,  knowledge  of  the  time  and  position  of  initial  dispersion 
and  a complete  knowledge  of  the  spatial  distribution  of  particles  at 
a later  time  admit  a relatively  simple  solution  of  the  problem.  In 
practice,  however,  one  seldom  has  such  complete  knowledge  of  the 
system.  Usually,  only  a sample  of  the  ensemble  is  observed  and  the 
time  and  place  of  initial  dispersion  may  not  be  known.  This  occurs, 
for  example,  in  the  cases  of  fragmentation  of  artificial  earth  satel- 
lites and  the  formation  and  evolution  of  meteor  streams. 

Thus,  the  study  of  our  inverse  problem  must  be  based  on  properties 
which  can  be  estimated  from  a subset  of  the  ensemble.  The  theory 
developed  in  the  sequel  is  based  on  the  moments  of  particle  position 
relative  to  a body  in  the  orbit  of  the  parent.  A truncation  at  some 
order  is  inevitable,  and  the  present  theory  is  truncated  at  the  fourth 
order  moments.  In  statistical  terminology,  the  moments  of  third  order 
incorporate  3kewness  of  the  distribution  while  fourth  order  moments 
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incorporate  kurtosis. 

Our  mathematical  treatment  of  the  inverse  problem  will  proceed  in 
stages  from  the  ideal  situation  of  complete  knowledge  of  , to  ? 

and  place  of  breakup  to  cases  of  progressively  less  observational 
material.  We  shall  deal  only  with  linear  dynamical  systems  in  this 
section.  Further,  let  us  consider  only  conservative  systems  so  that 
,t  ) = « This  is  done  Just  for  convenience  and  there  is 

no  difficulty  in  modifying  the  theory  for  non-zero  P . 

Given  f ( \ ) and  Time  and  Place  of  Initial  Disintegration 
To  begin  the  study  of  the  inverse  problem,  let  us  examine  the 
case  in  which  one  knows  the  spatial  density  completely  as  well  as  the 
time  of  initial  dispersion,  an  equivalent  statement  is  that  the  time 
t is  known.  From  equation  (21),  the  spatial  density  is  related  to  the 
momentum  distribution  function  according  to 


Mt  Yu  t )f 


Mio.t ) ^ - Yio,t)  V'(ji) 


(21  bis) 


We  seek  an  expression  for  f ' . The  simple  form  of  the  relation  (21) 
allows  us  to  write  immediately 


= | .iit  V(  o,  t )j  ^[m' \o,t ) ( £ -*■  Yi  j t ' V ' \ o t ) ^ ^ ) j 


(103) 


As  an  example,  suppose  were  known  to  be  of  ellipsoidal  form 
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f {-{  (%-p-  Ac^-  ^ , 


(l  oh) 


where  A is  a symmetric,  positive  definite  matrix  and  | is  a constant. 
The  momentum  distribution  which  produced  this  particle  distribution 
would  be  the  ellipsoidal  distribution 


(105) 


Q.=  [f  *Y'j,tiV'(o,t)^-M^].  M'Vo.t)  A +Y(o.t)V'Vo,t)|<-  v\s  J 

Given  Moments  of  ? ( ) and  the  Time  and  Place  of  Initial 

Disintegration 

Now  we  relax  the  requirement  on  knowledge  of  f ( ^ ) and  assume 
knowledge  of  the  low  order  moments  (order  zero  through  four).  A 
solution  of  the  inverse  problem  may  be  obtained  in  this  case  by 
introducing  expansions  for  p ( ^ ) and  G(g)  in  three-dimensional  Hermite 
polynomials  and  by  using  equation  (103)  to  relate  the  known  moments 
to  the  unknown  coefficients  in  the  expansion  for  G( g ) . There  is,  in 
principle,  no  necessity  to  truncate  the  expansion  at  fourth  order. 

The  development  given  below  may  be  carried  out  tr  any  desired  order. 

In  practice,  however,  it  seems  desirable  to  truncate  the  theory  at 
fourth  order. 
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To  begin  the  analysis,  let  us  recall  some  basic  properties  of 
three-dimensional  Hermite  polynomials  (Grad,  19^9) . The  polynomials 
are  defined  by 


with 


uj' ' d” 


Here  D*  is  the  differential  operator  of  all  partial  derivatives  of 
order  and  ^ is  a tensor  of  order  n whose  components  are  poly- 

nomials of  degree  n.  . Explicit  expressions  for  the  polynomials 
through  fourth  order  are: 


X'*'  = i 


» * • » 

i 

< i ) 


X. 


*J 


’fcJ  - 2ij 


v !3> 

x *j*.  = V*.  ~ 1 fJk.  <-  *-jJ ik  * ) 


+ ( *>)  <5x1  + * Su6)k  ) 


, x x,  i N 

h,  * t * j ' 


where  S is  the  Kronecker  delta. 

If  f<*}  is  a function  satisfying  the  condition 

eo 

J P (*)  u'1  il  X.OO  . 


(106) 


then  ( may  be  represented  by  the  Hermite  expansion 


A * aa  r" 


r /*- 


(10T) 


where  the  subscripts  on  and  a-1"'  are  n-indices.  Any  function  of 
compact  support,  hence  any  spatial  density  function,  satisfies  con- 
dition (106)  and  may  be  represented  by  the  series  (107).  Note  that 
not  all  terms  in  the  series  (107)  are  unique  because  *■'*'  occurs  once 
for  each  permutation  of  A . The  Hermite  coefficients  for  equation 
(107)  are  evaluated  from 


• \ C r . t*>  . 

L v = ~T  J ^ ^ v^>  ^ 

V r\  I v X - — 


(108) 


To  apply  these  concepts  to  the  solution  of  the  inverse  problem 
we  introduce  the  variable 


2 * A (j-*.)  , 


where  A is  a symmetric,  positive  definite  matrix  and  ^ is  a constant. 
A and  are  unspecified  at  this  point.  Next,  represent  G^l  by  the 
series 


t*.  ) +■  • • • 


(109) 


From  equation  (21),  it  follows  that  f(|>  is  represented  by  the  series 


(110) 
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where  here,  and  in  the  sequal,  the  time  dependent  matrices  are 


evaluated  at  the  arguments  t a , t ) 


The  next  step  is  to  relate  the  moments  of 


assumed  known,  to 


the  unknown  coefficients  a. 


The  zeroth  moment 


The  first  order  moments  involve  the  center  of  mass  of  the 


ensemble  and  we  obtain 


The  second  order  moments  involve  the  covariances  of  position 
relative  to  the  mean.  In  dyadic  notation,  we  obtain 


We  now  observe  that  the  expansion  (109)  for  G( £ ) is  valid  for 
any  positive  definite,  symmetric  matrix  A and  arbitrary  * , . 


Furthermore,  it  is  not  possible  to  separate  A andpa  from  a-‘'x  and 
*■  ^ . Therefore,  one  is  free  to  choose  A and  £,  and  determine  v , 

U> 

1 «.  or  vice  versa.  We  shall  follow  the  latter  course  and  choose 


— Ol.  — O 


. This  simplifies  the  moment  equations  considerably 


without  any  loss  of  generality. 

Before  proceeding,  let  us  ease  the  notational  burden  by  introduc- 


- " %“<•*> 


B - l AW\ )"' 


With  these  simplifications,  equations  (112)  and  (113)  become 


■*>=  L - VV’^.j  , 


(llU) 


^ — .tut  . ( M A I'A  ) =a«t  SB 

M *J  ij 


(115) 


The  third  and  fourth  order  moments  are  found  to  be 


llK.  » [>»./«."]  ZL  b.  b b <xt,J 

1 : J 77"t  -r  )>  k.  t <•  j t 


(116) 


<i<£££>j)lV  = b,,  b^b^  +!/*:>')  J*’, (117) 


where  b . & 

ll  'J 


To  complete  the  solution,  we  must  solve  equations  ( llU )-( 117 ) 


for 


A ) 

‘ M*. 


and  a. 


<*) 

iik: 


Let 


r "" 


where  X is  an  n-index.  From  equations  (llU)  and  (115)  we  obtain 


and 


Av  = ( na  r 11 ' m ■) - 1 


(118) 


(119) 


Note  that  these  relations  are  presaged  by  equations  ( 10U ) and  (105). 

To  invert  the  third  and  fourth  order  relations,  introduce  the  quantities 


Aw  = • B*' 

/ V ' 1 


Then  we  obtain 


7f  C 


(120) 


and 


. = 4. » j \ ^ ^iU  j ^ ji  * S i ^ 

1 ) *'■ 


(121) 


It  is  clear  from  equations  (120)  and  (121 ) that  higher  order 
moment  equations  would  present  no  inversion  difficulties.  At  all 
orders,  only  S is  needed  to  execute  the  inversion. 

Assembling  the  above  results,  we  obtain  the  following  expression 


6 5 


for  >»  <•  (i  > ; 


&■(  £ ) = ^ V | <A)  { ] t ji  ^ 0<.r  Pj ' t Pm.*  ^ - Z^rsi  tj  k. 


where 


uj  = {.-j  <*-*.>• 

F.  = *<■>>  - Ys7'’  , 

and 

a*  = ( M r‘”  m )_1  . 


Comments  on  More  General  Situations 

The  case  treated  in  the  previous  section  is  usually  appropriate 
for  the  analysis  of  satellite  breakups.  The  reason  for  this  is  that 
the  most  common  set  of  data  available  on  a breakup  is  a complete  set 
of  orbital  elements  for  each  fragment  and  there  are  several  independent 
methods  available  for  determining  the  time  and  place  of  breakup  from 
such  data. 

In  general,  if  one  knows  neither  the  time  nor  place  of  breakup, 
the  inverse  problem  can  be  solved  if  the  moments  of  f>  are  known  at 
two  different  times.  The  plausibility  of  this  statement  follows  from 
the  observation  that  Lambert's  theorem  provides  a method  for  calculating 
an  orbit  given  two  position  measurements  and  the  time  interval 
separating  the  measurements.  A method  of  solution  in  this  case  will  be 
outlined  below.  But  first  let  us  consider  a weaker  generalization. 

Suppose  we  have  two  sets  of  moments  associated  with  known  times. 
Suppose  we  know  neither  the  time  nor  place  of  breakup  but  do  know 
the  orbit  of  the  parent  body  from  which  the  fragments  originated.  This 
last  piece  of  information  relates  place  of  breakup  to  time  of  breakup. 
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applying  the  algorithm  developed  in  the  previous  section  to  both  sets 
of  moments  one  would  obtain  two  momentum  distribution  functions, 

and  <3\<  which  depend  on  the  unknown  time  of  breakup. 

Now  form  a "loss  function" 


<L  <t)  =•  J [&,<  fc  . (I22) 

Assuming  a unique  solution  to  the  inverse  problem,  one  may  find  the 
unknown  t by  finding  the  zero  of  1(0  . In  practice;  because  of 
observational  errors,  sampling  errors,  and  truncation  errors;  one 
should  seek  the  minimum  of  1.1*  > rather  than  a zero.  Once  t is  found, 
and  G0  will  agree  (or  nearly  so)  and  will  provide  the  solution. 

The  more  general  case  in  which  we  do  not  assume  knowledge  of  the 
parent  orbit  is  not  conceptionally  different.  In  this  case  the  two 
sets  of  moments  provide  two  momentum  distribution  functions  ' t 1 l » ) 
and  which  depend  on  the  unknown  time,  t,  and  the  unknown 

position  . The  appropriate  "loss  function"  is 


= J [<V  1 ip  • (123) 

Now,  minimizing  JL  with  repsect  to  the  four  parameters  ^-l‘t  will  yield 
the  solution. 

Both  of  the  general  cases  outlined  here  would  involve  considerable 
numerical  effort.  It  would  not  be  too  difficult  to  write  analytical 
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expressions  for  the  integrals  (122)  and  (123),  but  that  will  not  be 
undertaken  here.  Even  with  analytic  expressions  for  the  integrals, 
the  minimization  would  involve  a non-trivial  amount  of  computation. 
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